Для аналізу диференційної експресії генів було обрано клітини лінії А549 як біологічний об’єкт. Клітини лінії A549 — це клітини недрібноклітинної (NSCLC) аденокарциноми легень.
Дані було отримано з бази даних GEO (Gene Expression Omnibus). Завантажені дані описують експеримент з секвенуванням РНК (RNA-sequencing) на опромінених рентгенівським випромінюванням і контрольних клітинах A549.
Клітини A549 були піддані одноразовій дозі рентгенівського опромінення 2Гр або 4Гр за допомогою приладу X-RAD 160-225. Після опромінення клітини A549 (з одноразовою дозою опромінення), а також неопромінені клітини A549 були інкубовані ще 48 годин. Після цього було зібрано загальну РНК для високопродуктивного секвенування. Загальний дизайн експерименту: дослідження 3 умов по 3 репліки. Три групи: група 0Гр (контроль), група 2Гр та група 4Гр.
Відповідна стаття: doi: 10.1038/s41419-022-04561-x, PMID: 35190532.
kable(rse_colData_f |>
select(external_id, sra.study_title, sra.sample_title) |> head(), format = "html", table.attr = "class='table table-striped'")
| external_id | sra.study_title | sra.sample_title | |
|---|---|---|---|
| SRR8371688 | SRR8371688 | RNA-sequencing in irradiated and normal A549 cells. | 0GY-1 |
| SRR8371689 | SRR8371689 | RNA-sequencing in irradiated and normal A549 cells. | 0GY-2 |
| SRR8371690 | SRR8371690 | RNA-sequencing in irradiated and normal A549 cells. | 0GY-3 |
| SRR8371691 | SRR8371691 | RNA-sequencing in irradiated and normal A549 cells. | 2GY-1 |
| SRR8371692 | SRR8371692 | RNA-sequencing in irradiated and normal A549 cells. | 2GY-2 |
| SRR8371693 | SRR8371693 | RNA-sequencing in irradiated and normal A549 cells. | 2GY-3 |
kable(head(rse_rowData_f), format = "html", table.attr = "class='table table-striped'")
| source | type | bp_length | phase | gene_id | gene_type | gene_name | level | havana_gene | tag | |
|---|---|---|---|---|---|---|---|---|---|---|
| ENSG00000278704.1 | ENSEMBL | gene | 2237 | NA | ENSG00000278704.1 | protein_coding | BX004987.1 | 3 | NA | NA |
| ENSG00000277400.1 | ENSEMBL | gene | 2179 | NA | ENSG00000277400.1 | protein_coding | AC145212.2 | 3 | NA | NA |
| ENSG00000274847.1 | ENSEMBL | gene | 1599 | NA | ENSG00000274847.1 | protein_coding | AC145212.1 | 3 | NA | NA |
| ENSG00000276256.1 | ENSEMBL | gene | 2195 | NA | ENSG00000276256.1 | protein_coding | AC011043.1 | 3 | NA | NA |
| ENSG00000278198.1 | ENSEMBL | gene | 1468 | NA | ENSG00000278198.1 | protein_coding | AC011043.2 | 3 | NA | NA |
| ENSG00000273496.1 | ENSEMBL | gene | 1468 | NA | ENSG00000273496.1 | protein_coding | AC011841.1 | 3 | NA | NA |
kable(rse_assay_f |>
arrange(desc(SRR8371688)) |> head(), format = "html", table.attr = "class='table table-striped'")
| SRR8371688 | SRR8371689 | SRR8371690 | SRR8371691 | SRR8371692 | SRR8371693 | SRR8371694 | SRR8371695 | SRR8371696 | |
|---|---|---|---|---|---|---|---|---|---|
| ENSG00000198804.2 | 11438923 | 12624134 | 9876942 | 11286817 | 11500584 | 7659250 | 5559007 | 6564597 | 9430680 |
| ENSG00000115414.18 | 6535629 | 7642890 | 5454817 | 6763706 | 5559684 | 4748784 | 3533072 | 3375268 | 4997642 |
| ENSG00000198886.2 | 6206643 | 7357551 | 5660816 | 5944139 | 5889170 | 3652051 | 2882935 | 3346084 | 4630646 |
| ENSG00000165092.12 | 5324460 | 5907729 | 4752635 | 6099076 | 6003428 | 4324184 | 2738698 | 2820161 | 4443411 |
| ENSG00000198712.1 | 3950989 | 4529637 | 3600388 | 4185883 | 4218314 | 2704965 | 2025084 | 2605506 | 3313521 |
| ENSG00000198786.2 | 3844190 | 4525495 | 3422576 | 3768348 | 3622434 | 2426284 | 1887319 | 2046555 | 3086360 |
Побудовано графік щільності для окремого зразка (SRR8371688) та гістограму частоти значень експресії
# Створення таблиці метаданих
colData_new <- data.frame(
row.names = colnames(rse_filtered),
condition = factor(c(rep("0Gy", 3), rep("2Gy", 3), rep("4Gy", 3)))
)
# Матриця counts
countData <- assay(rse_filtered)
# Побудова графіку щільності
p1 <- ggplot(countData, aes(x = SRR8371688 )) +
geom_density(color = "darkblue", fill = "lightblue") +
xlim(0, 10e4)
count_summary <- countData |>
as.vector() |>
data.frame(counts = _) |>
group_by(counts) |>
summarise(frequency = n(), .groups = 'drop') |>
arrange(counts)
# Створення гістограми частоти значень
p2 <- ggplot(count_summary[1:200,], aes(x = counts, y = frequency, fill = frequency)) +
geom_bar(stat = "identity", width = 0.7) +
scale_y_log10(labels = scales::comma) +
labs(title = "Histogram of Counts vs Frequency (Logarithmic Scale)",
x = "Counts",
y = "Frequency (log10)") +
scale_fill_gradient(low = "#AA05BA", high = "#56B1F7") +
theme_minimal()
#p1
#p2
counts_info <- p1 + p2
ggsave("./figs/counts_info.png", plot = counts_info, width = 16, height = 6, dpi = 300)
knitr::include_graphics("./figs/counts_info.png")
Для аналізу диференційної експресії генів було застосовано DESeq Wald
Додатково було проведено коригування на хибний позитивний результат p.adj (FDR). Нижче наведено відповідний графік
res_wald_2Gy_vs_0Gy <- results(dds_wald, contrast = c("condition", "2Gy", "0Gy"))
res_wald_4Gy_vs_0Gy <- results(dds_wald, contrast = c("condition", "4Gy", "0Gy"))
res_wald_4Gy_vs_2Gy <- results(dds_wald, contrast = c("condition", "4Gy", "2Gy"))
hist_data_before <- bind_rows(
data.frame(padj = res_wald_2Gy_vs_0Gy$padj, Comparison = "2Gy vs 0Gy"),
data.frame(padj = res_wald_4Gy_vs_0Gy$padj, Comparison = "4Gy vs 0Gy"),
data.frame(padj = res_wald_4Gy_vs_2Gy$padj, Comparison = "4Gy vs 2Gy")
)
res_wald_2Gy_vs_0Gy$padj <- p.adjust(res_wald_2Gy_vs_0Gy$pvalue, method = "fdr")
res_wald_4Gy_vs_0Gy$padj <- p.adjust(res_wald_4Gy_vs_0Gy$pvalue, method = "fdr")
res_wald_4Gy_vs_2Gy$padj <- p.adjust(res_wald_4Gy_vs_2Gy$pvalue, method = "fdr")
hist_data_after <- bind_rows(
data.frame(padj = res_wald_2Gy_vs_0Gy$padj, Comparison = "2Gy vs 0Gy"),
data.frame(padj = res_wald_4Gy_vs_0Gy$padj, Comparison = "4Gy vs 0Gy"),
data.frame(padj = res_wald_4Gy_vs_2Gy$padj, Comparison = "4Gy vs 2Gy")
)
deseq_Wald_padj_b <- ggplot(hist_data_before, aes(x = padj, fill = Comparison, color = Comparison)) +
geom_density(alpha = 0.2, adjust = 1.5) + # Криві густини
scale_fill_manual(values = c("2Gy vs 0Gy" = "red", "4Gy vs 0Gy" = "blue", "4Gy vs 2Gy" = "forestgreen")) +
scale_color_manual(values = c("2Gy vs 0Gy" = "red", "4Gy vs 0Gy" = "blue", "4Gy vs 2Gy" = "forestgreen")) +
labs(title = "Розподіл коригованих p-значень", x = "p-значення", y = "Щільність") +
geom_vline(xintercept = 0.05, linetype = "dashed", color = "black") +
annotate("text", x = 0.05, y = Inf, label = "0.05", vjust = 1.5, hjust = -0.1, color = "black") +
theme_minimal()
deseq_Wald_padj_a <- ggplot(hist_data_after, aes(x = padj, fill = Comparison, color = Comparison)) +
geom_density(alpha = 0.2, adjust = 1.5) + # Криві густини
scale_fill_manual(values = c("2Gy vs 0Gy" = "red", "4Gy vs 0Gy" = "blue", "4Gy vs 2Gy" = "forestgreen")) +
scale_color_manual(values = c("2Gy vs 0Gy" = "red", "4Gy vs 0Gy" = "blue", "4Gy vs 2Gy" = "forestgreen")) +
labs(title = "Розподіл коригованих p-значень", x = "p-значення", y = "Щільність") +
geom_vline(xintercept = 0.05, linetype = "dashed", color = "black") +
annotate("text", x = 0.05, y = Inf, label = "0.05", vjust = 1.5, hjust = -0.1, color = "black") +
theme_minimal()
deseq_Wald_padj <- deseq_Wald_padj_b + deseq_Wald_padj_a
ggsave("./figs/deseq_Wald_padj.png", plot = deseq_Wald_padj, width = 16, height = 6, dpi = 300)
knitr::include_graphics("./figs/deseq_Wald_padj.png")
Для візуалізації у зміні експресії генів для випадку опромінення рентгеном 2Гр (порівняння з контролнм) та 4Гр (порівняння з контролнм) побудовано volcano plot. Окремо проводилось порівняння експресії генів при опроміненні 2Гр та 4Гр.
# png("./figs/deseq_Wald.png", width = 2400, height = 800)
# par(mfrow = c(1, 3), mar = c(5, 5, 4, 2))
# plotMA(res_wald_2Gy_vs_0Gy, ylim=c(-30,30), main="MA-plot: 0Gy vs 2Gy")
# plotMA(res_wald_4Gy_vs_0Gy, ylim=c(-30,30), main="MA-plot: 0Gy vs 4Gy")
# plotMA(res_wald_4Gy_vs_2Gy, ylim=c(-30,30), main="MA-plot: 2Gy vs 4Gy")
# dev.off()
# knitr::include_graphics("./figs/deseq_Wald.png")
res_wald_df_2Gy_vs_0Gy <- as.data.frame(res_wald_2Gy_vs_0Gy)
res_wald_df_2Gy_vs_0Gy$significant <- res_wald_df_2Gy_vs_0Gy$padj < 0.05
res_wald_df_4Gy_vs_0Gy <- as.data.frame(res_wald_4Gy_vs_0Gy)
res_wald_df_4Gy_vs_0Gy$significant <- res_wald_df_4Gy_vs_0Gy$padj < 0.05
res_wald_df_4Gy_vs_2Gy <- as.data.frame(res_wald_4Gy_vs_2Gy)
res_wald_df_4Gy_vs_2Gy$significant <- res_wald_df_4Gy_vs_2Gy$padj < 0.05
res_wald_df_2Gy_vs_0Gy$comparison <- "2Gy vs 0Gy"
res_wald_df_4Gy_vs_0Gy$comparison <- "4Gy vs 0Gy"
res_wald_df_4Gy_vs_2Gy$comparison <- "4Gy vs 2Gy"
res_wald_df_2Gy_vs_0Gy <- as.data.frame(res_wald_df_2Gy_vs_0Gy) |>
mutate(gene_id = rownames(res_wald_df_2Gy_vs_0Gy), comparison = "2Gy vs 0Gy") |>
left_join(gene_annotations, by = "gene_id")
# Оновлення res_df_4Gy_vs_0Gy
res_wald_df_4Gy_vs_0Gy <- as.data.frame(res_wald_df_4Gy_vs_0Gy) %>%
mutate(gene_id = rownames(.), comparison = "4Gy vs 0Gy") |>
left_join(gene_annotations, by = "gene_id")
# Оновлення res_df_4Gy_vs_2Gy
res_wald_df_4Gy_vs_2Gy <- as.data.frame(res_wald_df_4Gy_vs_2Gy) |>
mutate(gene_id = rownames(res_wald_df_4Gy_vs_2Gy), comparison = "4Gy vs 2Gy") |>
left_join(gene_annotations, by = "gene_id")
# Побудова volcano plot
plot_2Gy_vs_0Gy_W <- EnhancedVolcano(res_wald_df_2Gy_vs_0Gy,
lab = res_wald_df_2Gy_vs_0Gy$gene_name,
x = 'log2FoldChange',
y = 'padj',
xlim = c(-5, 5),
ylim = c(0, 10),
pCutoff = 0.05,
FCcutoff = 1.5,
pointSize = 5,
labSize = 6,
col = c('grey30', 'forestgreen', 'royalblue', 'red2'),
title = '2Gy vs 0Gy',
legendPosition = 'top')
plot_4Gy_vs_0Gy_W <- EnhancedVolcano(res_wald_df_4Gy_vs_0Gy,
lab = res_wald_df_4Gy_vs_0Gy$gene_name,
x = 'log2FoldChange',
y = 'padj',
xlim = c(-5, 5),
ylim = c(0, 10),
pCutoff = 0.05,
FCcutoff = 1.5,
pointSize = 5,
labSize = 6,
col = c('grey30', 'forestgreen', 'royalblue', 'red2'),
title = '4Gy vs 0Gy',
legendPosition = 'top')
plot_4Gy_vs_2Gy_W <- EnhancedVolcano(res_wald_df_4Gy_vs_2Gy,
lab = res_wald_df_4Gy_vs_2Gy$gene_name,
x = 'log2FoldChange',
y = 'padj',
xlim = c(-5, 5),
ylim = c(0, 10),
pCutoff = 0.05,
FCcutoff = 1.5,
pointSize = 5,
labSize = 6,
col = c('grey30', 'forestgreen', 'royalblue', 'red2'),
title = '4Gy vs 2Gy',
legendPosition = 'top')
volc_wald_plot <- plot_2Gy_vs_0Gy_W + plot_4Gy_vs_0Gy_W + plot_4Gy_vs_2Gy_W + plot_layout(nrow = 1)
#print(volc_wald_plot)
ggsave("./figs/volc_wald.png", plot = volc_wald_plot, width = 20, height = 6, dpi = 300)
ggsave("./figs/volc_plot_2Gy_vs_0Gy.png", plot = plot_2Gy_vs_0Gy_W, width = 20, height = 20, dpi = 300)
knitr::include_graphics("./figs/volc_plot_2Gy_vs_0Gy.png")
ggsave("./figs/volc_plot_4Gy_vs_0Gy.png", plot = plot_4Gy_vs_0Gy_W, width = 20, height = 20, dpi = 300)
knitr::include_graphics("./figs/volc_plot_4Gy_vs_0Gy.png")
ggsave("./figs/volc_plot_4Gy_vs_2Gy.png", plot = plot_4Gy_vs_2Gy_W, width = 20, height = 20, dpi = 300)
knitr::include_graphics("./figs/volc_plot_4Gy_vs_2Gy.png")
Для візуалізації перекриття між списками генів, які є статистично значущими (з коригованим p-значенням < 0.05) та мають великий зсув логарифмічних змін (log2FoldChange > 1.5 для підвищеної експресії та log2FoldChange < -1.5 для зниженої експресії) у різних порівняннях, використали upset plot.
res_2Gy_0Gy_up <- res_wald_df_2Gy_vs_0Gy |> as.data.frame() |> filter(pvalue < 0.05 & log2FoldChange > 1.5)
res_2Gy_0Gy_down <- res_wald_df_2Gy_vs_0Gy |> as.data.frame() |> filter(pvalue < 0.05 & log2FoldChange < -1.5)
res_4Gy_0Gy_up <- res_wald_df_4Gy_vs_0Gy |> as.data.frame() |> filter(pvalue < 0.05 & log2FoldChange > 1.5)
res_4Gy_0Gy_down <- res_wald_df_4Gy_vs_0Gy |> as.data.frame() |> filter(pvalue < 0.05 & log2FoldChange < -1.5)
res_4Gy_2Gy_up <- res_wald_df_4Gy_vs_2Gy |> as.data.frame() |> filter(pvalue < 0.05 & log2FoldChange > 1.5)
res_4Gy_2Gy_down <- res_wald_df_4Gy_vs_2Gy |> as.data.frame() |> filter(pvalue < 0.05 & log2FoldChange < -1.5)
venn_list <- list("2Gy vs 0Gy(up)" = res_2Gy_0Gy_up$gene_name,
"2Gy vs 0Gy(down)" = res_2Gy_0Gy_down$gene_name,
"4Gy vs 0Gy(up)" = res_4Gy_0Gy_up$gene_name,
"4Gy vs 0Gy(down)" = res_4Gy_0Gy_down$gene_name,
"4Gy vs 2Gy(up)" = res_4Gy_2Gy_up$gene_name,
"4Gy vs 2Gy(down)" = res_4Gy_2Gy_down$gene_name)
ven <- venndetail(venn_list)
png("./figs/upset_diagram.png", width = 1600, height = 1600, res = 300)
plot(ven, type = "upset")
dev.off()
Наступним кроком було проведено PCA аналіз. Scree plot Червона лінія відображає кумулятивну суму варіацій, яка показує, що 80% варіації даних припадає на п’ять головних компонент. При цьому PC1 пояснює близько 50% варіації. Для того, щоб визначити кількість головних компонент, які варто зберегти для аналізу були використані метод Горна та elbow, які зображені на графіку штриховою лінією. Biplot комбінує інформацію про зразки і змінні в одному графіку. По осях зображені дві головні компоненти: PC1 і PC2, кожна з яких пояснює відсоток варіації в даних: PC1 — 44.98%, PC2 — 12.94%. Точки відповідають різним спостереженням за умов випромінювання X-хвилями 0Гр, 2Гр, 4Гр, які приблизно кластеризуються відповідно умовам впливу. Щоб перевірити чи є “правильна” кластеризація експериментів, був проведений кластерний аналіз методом k-means з 5-ма головними компонентами, адже саме стільки компонент пояснюють більшу варіацію даних. Loadings plot показуює, наскільки кожен ген впливає на кожну головну компоненту. Це значення визначає, наскільки цей ген сприяє поясненню варіації в даних.
vst_data <- assay(vst(dds_wald))
pca_result <- PCAtools::pca(vst_data, metadata = colData_new, removeVar = 0.1)
horn <- parallelPCA(vst_data)
elbow <- findElbowPoint(pca_result$variance)
var_explained <- which(cumsum(pca_result$variance)> 80)[1]
pscree <- PCAtools::screeplot(pca_result,
axisLabSize = 18,
titleLabSize = 22,
vline = c(horn$n, elbow, var_explained), hline = 85,
returnPlot = FALSE) +
geom_label(aes(x = horn$n, y = 60,
label = 'Horn\'s', vjust = -1, size = 8)) +
geom_label(aes(x = elbow, y = 60,
label = 'Elbow method', vjust = -3, size = 8)) +
geom_label(aes(x = var_explained , y = 85,
label = '80% explained variation', vjust = -1, size = 8))
ppairs <- PCAtools::pairsplot(pca_result,
hline = 0, vline = 0,
pointSize = 2,
gridlines.major = FALSE, gridlines.minor = FALSE,
colby = 'condition',
plotaxes = TRUE,
margingaps = unit(c(0.01, 0.01, 0.01, 0.01), 'cm'),
returnPlot = FALSE)
# Створюємо нові назви для рядків, зіставляючи gene_id з gene_name
new_row_names <- gene_annotations$gene_name[match(rownames(pca_result$loadings), gene_annotations$gene_id)]
# Якщо є дублікат імен, додамо індекси для унікальності
new_row_names <- make.unique(new_row_names)
# Перезаписуємо назви рядків у pca_result$loadings
rownames(pca_result$loadings) <- new_row_names
pbiplot <- PCAtools::biplot(pca_result,
axisLabSize = 18,
# loadings parameters
showLoadings = TRUE,
sizeLoadingsNames = 4,
colLoadingsNames = 'red4',
#points
lab = NULL,
colby = 'condition',
colLegendTitle = 'Condition',
colkey = c('0Gy' = '#228B22', '2Gy' = '#7B68EE', '4Gy' = '#FF8C00'),
pointSize = 5,
#plotview
hline = 0,
vline = c(-5, 0, 5),
vlineType = c('dotdash', 'solid', 'dashed'),
gridlines.major = FALSE, gridlines.minor = FALSE,
legendPosition = 'right',
legendLabSize = 16,
legendIconSize = 8.0,,
subtitle = 'PC1 versus PC2',
caption = '5 PCs ≈ 80%',
returnPlot = FALSE)
ploadings <- PCAtools::plotloadings(pca_result,
axisLabSize = 18,
subtitle = 'PC1, PC2, PC3, PC4, PC5',
labSize = 3,
returnPlot = FALSE)
pca_df <- as.data.frame(pca_result$rotated)
# Додаємо інформацію про condition та інші метадані
pca_df$condition <- colData_new$condition
# Використовуємо більше компонент для кластеризації, наприклад, перші 5
set.seed(1234)
km <- kmeans(x = pca_df[, 1:5], centers = 3, nstart = 25) # Використовуємо 5 компонент і 25 запусків
# Додаємо результати кластеризації в датафрейм
pca_df$cluster <- as.factor(km$cluster)
# Візуалізація
pkmean <- autoplot(km, data = pca_df, frame = TRUE, frame.type = 'norm', size = 4, alpha = 0.6) +
aes(shape = condition, color = cluster) +
scale_shape_manual(values = c('0Gy' = 16, '2Gy' = 17, '4Gy' = 18)) +
labs(subtitle = 'with 5 PCs',
color = 'cluster', shape = 'Condition') +
coord_equal() +
theme_minimal(base_size = 18)+
theme(legend.position = 'right', # Позиція легенди
panel.border = element_rect(color = 'black', fill = NA, size = 1)) # Додання рамки
# Додаємо заголовки до кожного графіка
pscree <- pscree + ggtitle("A Scree plot") + theme(plot.title = element_text(face = "bold", hjust = -0.01, size = 18))
ploadings <- ploadings + ggtitle("B Loading plot") + theme(plot.title = element_text(face = "bold", hjust = -0.01, size = 18))
pbiplot <- pbiplot + ggtitle("C Bi-plot") + theme(plot.title = element_text(face = "bold", hjust = -0.01, size = 18))
pkmean <- pkmean + ggtitle("D K-means Clusters") + theme(plot.title = element_text(face = "bold", hjust = -0.01, size = 18))
ppairs <- ppairs + ggtitle("E Pairs plot") + theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 18))
# Комбінування графіків у сітку
combined_plot <- (pscree | ploadings) / (pbiplot | pkmean) / ppairs
ggsave("./figs/combined_plot.png", combined_plot, width = 20, height = 30)
knitr::include_graphics("./figs/combined_plot.png")
Після проведення аналізу диференційної експресії було визначено роль цих генів в клітинах. Для цього було проведено аналіз збагачення за допомогою Enrichr
knitr::include_graphics("./figs/enrichr_results_all_up.png")
knitr::include_graphics("./figs/enrichr_results_2Gy_0Gy_up.png")
knitr::include_graphics("./figs/enrichr_results_4Gy_0Gy_up.png")
knitr::include_graphics("./figs/enrichr_results_4Gy_2Gy_up.png")
Для того, щоб явно побачити роль генів з підвищеною та зі зниженою експресією у клітинах лінії А549 було проведено аналіз метаболічних шляхів циклу клітини, проліферації, апоптозу.
logFC <- path_4Gy_vs_0Gy$log2FoldChange
names(logFC) <- path_4Gy_vs_0Gy$entrezgene_id
pathview(gene.data = logFC,
pathway.id = "hsa05207",
species = "hsa",
#kegg.native = T
ut.suffix = "Non-small cell lung cancer (специфічно для A549)")
pathview(gene.data = logFC,
pathway.id = "hsa04110",
species = "hsa",
#kegg.native = T
ut.suffix = "Cell cycle")
pathview(gene.data = logFC,
pathway.id = "hsa04115",
species = "hsa",
#kegg.native = T
ut.suffix = "p53 signaling pathway")
pathview(gene.data = logFC,
pathway.id = "hsa04725",
species = "hsa",
#kegg.native = T
ut.suffix = "NAChR")
pathview(gene.data = logFC,
pathway.id = "hsa04210",
species = "hsa",
#kegg.native = T
ut.suffix = "Apoptosis")
knitr::include_graphics("hsa05207.pathview.png")
knitr::include_graphics("hsa04110.pathview.png")
knitr::include_graphics("hsa04115.pathview.png")
knitr::include_graphics("hsa04210.pathview.png")
knitr::include_graphics("hsa04725.pathview.png")
## {-}